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We discuss results from 3+1-D numerical simulations of SU(2) Yang-Mills equations for 
an unstable Glasma expanding into the vacuum after a high energy heavy ion collision. We 
expand on our earlier work on a non-Abelian Weibel instability in such a system and study 
the behavior of the instability in greater detail on significantly larger lattices than previously. 
We establish the time scale for the onset of the instability and demonstrate that the growth 
■ rate is robust as one approaches the continuum limit. For large violations of boost invariance, 

{N| \ non-Abelian effects cause the growth of soft modes to saturate. At late times, we observe 

significant creation of longitudinal pressure and a systematic trend towards isotropy. These 
^ \ time scales however are significantly larger than those required for early thermalization in 

heavy ion collisions. We discuss additional effects in the produced Glasma that may speed 
up thermalization. 



O 



o 
o 



I. INTRODUCTION 

An outstanding theoretical puzzle in high energy heavy ion collisions is to demonstrate the ther- 
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malization of the quark-gluon matter produced in these collisions. Experiments at the Relativistic 
Heavy Ion Collider (RHIC) at Brookhaven National Laboratory (BNL) indicate that such a ther- 
malized state, the Quark Gluon Plasma (QGP), has been formed in collisions of ultrarelativistic 
nuclear beams 0] . Understanding thermalization from first principles in Quantum Chromodynam- 
ics (QCD) is complicated by the interplay of several overlapping time scales. The collision itself, 
Q_i| in the framework of the Color Glass Condensate (CGC) effective theory 0], can be understood as 

the collision of coherent classical Yang-Mills fields 0, . The typical momenta of these fields are 
characterized by the scale 1 Q s 0,0]. Because it is the only relevant time scale in the problem, the 
formation time of gluons after the collision is of order 1/Q S 0,E1- The nuclei are highly Lorentz 
contracted - the pressure gradients in the longitudinal direction after the collision are therefore 
^ enormous and the system expands outwards at nearly the speed of light. The initial space-time evo- 

lution of the produced classical fields is described by solutions of Yang-Mills equations with CGC 
initial conditions 0,0,0,0]. The produced gluons begin to scatter; the strength of their scattering 
is determined by an infrared Debye mass ran which screens the range of their interactions. The 
occupation number / of produced gluons, initially of order of the inverse strong coupling, / ~ l/a s , 
decreases with time. It was believed initially that 2 — > 2 ^lj and later 2 — > 3 [pa Ipi ] scattering 
processes could thermalize the system. Parametrically, however, these scattering processes take 
a long time of order r t herm. ~ 13/5 rr i n the "bottom up" thermalization scenario 13]. Matter 
in this pre-equilibrium phase, where several time scales compete, has been termed a "Glasma" 
to describe its transitory behavior from coherent Color Glass fields to thermalized Quark Gluon 
Plasma Q. 

It has been suggested for some time that instabilities, analogous to the Weibel instability [1_ 
in plasma physics, may play an important role in thermalization of the Glasma Recently, 
the specific mechanism in the "Color Glass/Bottom Up" scenario triggering the instability was 



1 In the CGC framework, this scale is simply related, in leading order, to /j 2 , the color charge squared per unit area: 
Ql = g 4 fJ 2 N c \n(g 2 u/A)/2n, where A is an infrared cut-off. 
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identified 13] as arising from the change in sign of the Debye mass squared for anisotropic mo- 
mentum distributions [la]. One can view this, in the configuration space of the relevant fields, as 
the development of specific modes for which the effective potential is unbound from below [l9(. 
Detailed simulations in the Hard-Loop effective theory in 1+1-dimensions 12, HI 3 and in 3+1- 
dimensions 21 , z| have confirmed the existence of this non-Abelian Weibel instability. Particle- 
field simulations of the effects of the instability on thermalization have also been performed re- 
cently [23J, [24|. 

All these simulations consider the effect of instabilities in systems at rest. However, as discussed 
previously, the Glasma expands into the vacuum at nearly the speed of light. Recently, we pre- 
sented first results on 3+1-D numerical simulations of Yang-Mills equations which demonstrated 
the existence of a non-Abelian Weibel instability in the expanding Glasma [2j|. Such an instability 
was not seen in previous numerical simulations of the Glasma E3] which assumed strict boost 
invariance, and therefore, dynamics in 2+1 dimensions. Remarkably, arbitrarily weak violations of 
boost invariance trigger the non-Abelian Weibel instability. Another striking feature of our simu- 
lations was that the unstable fields in the expanding Glasma grow proportional to exp(Ty/ g 2 fir) as 
opposed to the usual exponential form. The former functional form was predicted by the authors of 
Ref . 0] . They deduced it simply from the fact that the scale for the growth rate in the expanding 
system is set by the Debye mass which depends on the proper time r(= yt 2 — z 2 ) as mo oc 1/y/r. 

In this work, we present more details and expand on the results presented in Ref. j2f|. In 
particular, we look at much larger transverse and longitudinal lattices and study the dependence of 
the growth rate on the volume and continuum limits. Further, we extend the studies of Ref. jj^J to 
much larger violations of boost invariance than considered previously. While kinematic violations 
of boost invariance are small, dynamical small x quantum evolution effects contribute to much 
larger violations of boost invariance. We consider one particular realization of such violations of 
boost invariance. We are able to follow all the stages of the evolution: we determine the time scale 
when the exp(T\/<? 2 //T) growth of unstable modes starts to dominate over the initial spectrum, the 
saturation of the growth of these modes and the subsequent generation of longitudinal pressure as 
the system evolves towards isotropy. The time scales we obtain for this last stage are much larger 
than natural time scales for heavy ion collisions. 

Recently, a field theory formalism was developed to compute multiplicities in field theories with 
strong time dependent sources j (with j ~ 1/g, where g is the QCD coupling). The prototype 
for such a theory is the CGC jig]. In this framework, the results of Refs. H, 0, 0] can be 
understood as a leading order contribution to the inclusive multiplicities, while the computation 
performed here can be seen as an (approximate) means of computing a piece of the next-to-leading 
order multiplicities. Other contributions are not included in this computation. (We will elaborate 
on these remarks in section V.B.) Alternately, the results of Ref. 26] can be formulated as a 



Boltzmann equation with a source term 23] - the full next-to-leading order contributions will 



provide significant contributions towards thermalization. Our computations here must therefore 
be viewed as an essential but not exclusive piece of a complete computation of the evolution of the 
Glasma into the QGP. 

This paper is organized as follows. In section II, we discuss the formulation of the problem 
of nuclear collisions in the CGC framework. In particular, we discuss how violations of boost 
invariance are implemented. We discuss details of the lattice formulation of the problem in the 
following section. In section IV, we discuss numerical results for the non-Abelian Weibel instability. 
In particular, we study how the Fourier modes of the longitudinal pressure grow as a function of 
proper time. We show explicitly the distribution of unstable longitudinal (k||k z ) modes and their 
evolution. We note a very interesting behavior for the hardest momentum mode that is still 
unstable, is max : while otherwise growing slowly, f max increases dramatically when the maximum 
amplitude of an unstable mode reaches a certain critical value. A plausible interpretation of this 
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effect is that when the amplitude of the longitudinal fluctuations - which roughly translates as 
the amplitude of the unstable transverse magnetic field modes - become large enough, the Lorentz 
force on transverse (k||kj_) "particle" modes of the gluon fields is sufficient for them to acquire 
significant longitudinal momenta. Hence the rapid increase in v max . This effect is enhanced for 
large boost invariance violating "seeds", which are discussed further in section V. Because this 
"Lorentz force effect" is significant for large seeds, we notice a rise in the longitudinal pressure, 
albeit the effect becomes measurable only at late times. The increase in the longitudinal pressure 
is accompanied by a decrease in the transverse pressure; this supports our conjecture that hard 
transverse gluon modes display a significant change in their trajectories when the modes of the 
transverse magnetic fields are large enough. We observe a clear trend towards isotropy and quantify 
the change in the dependence of the energy density with proper time. Despite this trend, the time 
scales over which these effects occur are much too large to explain the early thermalization at 
RHIC. This may be because, as discussed previously, our numerical simulations may not include 
all the processes necessa ry f or thermalization. In the final section, we summarize our results and 
discuss work in progress 28] on thermalization in the Glasma. 

II. NUCLEAR COLLISIONS IN THE COLOR GLASS CONDENSATE 

We will provide here a brief review of nuclear collisions in the Color Glass Condensate framework. 
A more detailed discussion can be found in Ref. 0. We will first discuss the 2+1-dimensional boost 
invariant formulation of the problem before extending our discussion to the more general 3+1-D 
problem. 

A. Gluon production from classical fields 

In nuclear collisions at very high energies, the hard valence parton modes in each of the nuclei 
act as highly Lorentz contracted, static sources of color charge for the soft wee parton, Weizsacker- 
Williams modes in the nuclei. By hard and soft, we mean large x or small x, where x is the 
longitudinal momentum fraction of partons in the colliding nuclei. Soft x modes can be understood 
as modes that are coherent across the longitudinal extent of the nucleus, or equivalently, x <C A^ 1 / 3 . 
With increasing energy, the scale separating soft and hard modes shifts towards smaller values of x. 
How the sources are modified with this changing scale is quantified by a Wilsonian RG procedure-a 
discussion and relevant references can be found in Ref. |2j. In a nuclear collision, the hard sources 
are described by the current 

j^a = fii* ft(x ± )8(x-) + p%(x ± )8(x+) , (1) 

where the color charge densities pf 2 of the two nuclei are independent sources of color charge 
on the light cone. Note that x^ = (t ± z)/2. The 8 functions represent the fact that Lorentz 
contraction has squeezed the nuclei to infinitesimally thin sheets. The absence of a longitudinal 
size scale ensures that the gauge fields generated by these currents will be boost-invariant, namely, 
independent of the space time rapidity rj = atanhf. The gauge fields before the collision are 
obtained by solving the Yang-Mills equations D^F^ U = j", with 

= 8^ + ig[A^ .} , = d^Ay - 8 U A^ + ig[Ap,A v ] , (2) 

the gauge covariant derivative and field strength tensor, respectively, in the fundamental represen- 
tation, and [Ap, .] denotes a commutator. 



4 



Gluon distributions are simply related to the Fourier transform A p (k±) of the solution of the 
Yang-Mills equations by {A fl (k±)A fl (k±)) p . The averaging over the classical charge distributions 
is defined by the expression 

(0) p = J d Pldp2 0( Pl , P2 ) exp (-/dV ^^^ ) , (3) 

and is performed independently for each nucleus with equal Gaussian weight P 2 . Here P 2 oc A 1 / 3 
fm~ 2 denotes the color charge squared per unit area in each of the (identical) nuclei with atomic 
number A. For very large nuclei, P is much larger than the fundamental QCD scale, fj? » Aq CD 
- this asymptotic condition justifies our applying weak coupling techniques to the problem. Such 
a Gaussian weight is justified [3. liol. EH in the limit of A 3> 1 and a s Y <C 1, providing a window 
ln(^4 1 / 3 ) < Y < A 1 /^ in rapidity 2 Y, which exists for large A. This window of applicability, for 
large nuclei, can be extended to larger values of Y by using the Balitsky-Kovchegov equation 
to evolve the sources px,2 to higher rapidities. This small x renormalization- group evolution of 
the source densities preserves the Gaussian structure of the sources - however, it is now non-local 
and /i 2 — > /i 2 (x± — y± ) . We do not expect this generalization to qualitatively modify the results 
discussed here. We will return to this discussion when we discuss violations of boost invariance. 
Before the nuclei collide (t < 0) , a solution of the equations of motion is 0] 

A ± = ; A* = 9 e {x-)9 e (-x + )a\(x±) + 9 e (x + )8 e {-x~)ai{x ± ) , (4) 

where, here and in the following, the transverse coordinates x,y have been collectively labeled by 
the Latin index i. The e subscripts on the ^-functions denote that they are smeared by an amount 
e in the respective x^ light cone directions. We require that the functions a l m (x±) (m = 1, 2 denote 
the labels of the colliding nuclei) satisfy F %3 = 0, namely, that they be pure gauge solutions of the 
equations of motion. The functions a l m satisfy 

- D i a \m) =P(m)(xl). ( 5 ) 
This equation has an analytical solution given by 0, 

« l (m) = — e iA <™> d'e-^ , Vi A (m) = -g p [m) . (6) 

To obtain this result one has to assume path ordering in x ± respectively for nuclei 1 and 2; we 
assume that the limit e — > is taken at the final step. 

We now introduce the proper time r = \/t 2 — z 2 = \/2x + x~ - the initial conditions for the 
evolution of the gauge field in the collision are formulated on the proper time surface r = 0. They 
are obtained by generalizing the previous ansatz for the gauge field to 

A i {x-,x + ,x ± ) = Q £ (x~)Q £ (-x + )a\(x ± ) + e e (-x~)e € (x + )al(x ± ) 
+@ e (x~)@ e (x + )ai(x~ , x + , x L ) 
A ± = ±x ± e e (x-)G € (x + )(3(x- ,x+ ,x ± ), (7) 

where we adopt the Fock-Schwinger gauge condition A T = x + A~ + x~ A + = 0. This gauge is 
an interpolation between the two light cone gauges A^ on the x ± = surfaces respectively. The 
sole purpose of the ansatz for the gauge fields here is to obtain the unknown gauge fields a^, /3 



Here we refer to the momentum space rapidity Y = i ln( ) = Ybcam — m (^ 
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in the forward light cone from the known gauge fields ckm 2) of the respective nuclei before the 
collision. This is achieved by invoking a physical "matching condition" which requires that the 
Yang-Mills equations D^F^ V = J v be regular at r = 0. The <5-functions of the current in the 
Yang-Mills equations therefore have to cancel identical terms in spatial derivatives of the field 
strengths. Interestingly, it leads to the unique solution 0] 

a\(x + ,x~,x ± ) = a\(x L ) + a 2 (x ± ), P(x + ,x~,x ± ) = --ig [ai t i(x±), a 2 (x_L)]. (8) 

Further, the only condition on the derivatives of the fields that would lead to regular solutions 
are d r a\ r= o , d T a^_| T= o = 0. This completes the derivation of the initial conditions in the boost- 
invariant case. 

Gyulassy and McLerran 33] argue that the equations of motion with the above boundary 
conditions remain unchanged even when the fields a\ 2 before the collision are smeared out in 
rapidity to properly account for singular contact terms in the equations of motion. We concur and 
will discuss violations of boost invariance later. 



B. Hamiltonian Chromodynamics in r, r\ coordinates 



While the Yang-Mills equations discussed above can be solved perturbatively 0, 0] i n the 
limit a s fi <C k±, it is unlikely that a simple analytical solution exists in general. The classical 
solutions have to be determined numerically for r > 0. The straightforward procedure would be to 
discretize the Yang-Mills equations. However, gauge invariance is best ensured by solving Hamil- 
ton's equations on the lattice. In practice therefore, we need to construct the lattice Hamiltonian 
and obtain the corresponding lattice equations of motion. 

The gluonic part of the QCD action in general coordinates takes the form 



S = — — J dTdt]dx± \J — detg^Tr 



dTdr]dx±j£, 



0) 



where for convenience we choose to keep the gauge-covariant but not coordinate-covariant F^v 
defined by Eq. ©. For the purpose of solving the Yang-Mills equations for a heavy-ion collision 
on a lattice, we shall work in the light cone coordinates r = V2x + x~ and r\ = \ ln(^). Our gauge 
condition x + A~ + x~A + = transforms to A T = 0, and the Lagrangian density becomes 



£ = rTr 



F 2 

rr\ 



F 2 



F 2 



3r\A-q 



(10) 



due to the non-trivial metric g^ u = diag(l, —1, —1, — r 2 ). Note that the contribution of the hard 
valence current to C drops out in the boost-invariant case 3 and can be neglected as long as one 
only considers a small region around central space-time rapidities r] = 0. The conjugate momenta 
are then 

dC „ . „ dC 1 



rd T Ai 



EL 



d(d T A, 

which we use to construct the Hamiltonian density 
H = Ei{d T Ai) + Er,(d T A v ) - C = Tr 



d{d T Ar,) 



d T A r . 



rp2 7?2 

— H — — + tE 2 + tF 2 

T T n xy 



(11) 



(12) 



3 As we shall see, the dependence on the source densities is contained in the initial conditions through the dependence 
of the gauge fields on the source densities. 
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Using finally 



on 

dE, 



dU 

dA,. 



-d T E, 



(13) 



one finds the equations of motion in the Hamiltonian formulation, 



drAi = ^, 
T 



d T A.ri 



The Gauss law constraint is simply 

DiEi + D V E V = . 



d T Erj 



TEr, 

t~ x DaF.. 



3 JV 



(14) 



(15) 



If the sources are strictly 5- function sources on the light cone, the Yang-Mills equations are 
independent of the space-time rapidity. One therefore has 



Ai(x ± ,r],T) = Ai(x±,r) ; A v = <^(x±,t) 



(16) 



which results in F 7 
case to 



■if i 



-Di$. The Hamiltonian in Eq. fi2j) then reduces in the boost invariant 



H = Ei{d T Ai) + E v (d T A v ) - C = TV 



El (B^f - „ 
-S- + V ' + rEl + tF 2 

T T 



xy 



(17) 



which is the QCD Hamiltonian in 2+1-dimensions coupled to an adjoint scalar field. Boost invari- 
ance clearly greatly simplifies the numerical simulations of these equations of motion. 

The initial conditions for the gauge-fields in the boost-invariant case from Eq. (jHJ become 

A l (x±) = am(ii) + afa(x±), A n (x±) = 0, Ei(x±) = 0, E v (x±) = ig [a l (1) , a ? (2) ]. 

(18) 

The magnetic fields being defined as Bk = ek^E^ , these initial conditions suggest that B, n / 
and Bi = 0. Note that the latter condition follows from the constraint on the derivatives of the 
gauge field that ensure regular solutions at r = 0. Thus one has E V ,B V ^ and large, as well 
as Ei,Bi = 0. This is in sharp contrast to the electric and magnetic fields of the nuclei before 
the collision - they are purely transverse! Some of the consequences of these initial conditions 
were discussed previously by Kharzeev, Krasnitz and Venugopalan [^Hl; their importance was 
emphasized recently by Lappi and McLerran 



C. Violations of Boost-Invariance 

In heavy-ion collisions, one clearly does not have exact boost-invariance. Besides simple geomet- 
ric effects (a heavy nucleus can never be Lorentz contracted into an infinitely thin sheet), one also 
has to take into account quantum fluctuations at high energies. These lead to violations of boost 
invariance that are of order unity over rapidity scales Y ~ ctj 1 . Consequently, a more realistic 
model of a heavy-ion collision will have color sources that have a finite width in x ± instead of the 
idealized 5 functions of Eq. (|T|). Nevertheless, the nuclei are highly localized on the light cone, and 
the same matching conditions required for regular solutions in the forward light cone apply. 

We will assume here that the initial conditions in Eq. (JSJ) can be extended to the boost non- 
invariant case. Ignoring details of exact matching of the Yang-Mills equations on the light cone, 
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the qualitative difference with respect to idealized boost-invariant case are rapidity fluctuations due 
to violations of boost invariance. In what follows, we will study two models of initial conditions 
containing rapidity fluctuations 4 . We construct these by modifying the boost-invariant initial 
conditions in Eq. (|18|) to 

Ei(x±,rj) = 8Ei(x±,r]), E rj (x±,rj) = ig[a i w ,a\ 2) ]+5E ri (x±,r]) (19) 

while keeping Ai,A v unchanged. The rapidity dependent perturbations dEi^SE^ are in principle 
arbitrary, except for the requirement that they satisfy the previously mentioned Gauss law. For 
these initial conditions, it takes the form 

Di5Ei + D ri E v = . (20) 

We construct these perturbations as follows: 

i) we first generate random configurations SEi(x±) with (5Ei(x±)5Ej(y±)) = 5ij5(x± — y±). 

ii) Next, for our first model of rapidity perturbations, we generate a Gaussian random function 
F{rj) with amplitude A 

< F{rj)F{r}') >= A 2 6(r) - rf). (21) 

For the second model, which we shall discuss further in section^ we generate the Gaussian random 
function but then remove the high-frequency components of F{rf) by applying a "band-filter" with 
strength b: 

F(rj) — > F(rj) = J — exp(— ivrj) exp(— \v\b) J drj exp(irjv)F(r}) . (22) 

The white noise Gaussian fluctuations of the previous model ensure that the amplitudes of all 
modes are of the same size. The high momentum modes dominate bulk observables such as the 
pressure. The instability we will discuss shortly is sensitive to the infrared modes but its effects are 
obscured by the higher momentum modes. This is particularly acute for large violations of boost 
invariance. Damping the high frequency modes of the white noise spectrum will therefore allow us 
to also study larger values of A, or "large seeds" that violate boost-invariance. 
hi) For both models, once F{rj) is generated, we then obtain 

5Ei{xi,rj) = d n F(n)5Ei{x x ), E v ( Xi ,rj) = -F^D^E^x^) . (23) 

This construction explicitly satisfies Gauss' law. 

We note that in order to implement rapidity fluctuations in the above manner, one has to have 
r > 0. This is a consequence of the chosen coordinates and does not have a physical origin. We 
will therefore will use these initial conditions at r = r m it with < -C Qj 1 and show our results 
are fairly independent on the specific choice of r m it. 



III. LATTICE SIMULATIONS 

The Hamiltonian formulation in the last section is very convenient since it allows for a straight- 
forward discretization that can be used to solve the Yang-Mills equations on a space-time lattice in 



4 These models of violations of boost invariance are only weakly motivated at present. The only condition imposed 
is that these satisfy Gauss' law. A more complete theory should specify, from first principles, the initial conditions 
in the boost non-invariant case. 
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(r,x,y,rj). We follow previous discussions of lattice equations of motion and initial conditions 
For the reader's convenience, these are reproduced in Appendix EI We work on a lattice that has 
periodic boundary conditions in the transverse plane as well as in space-time rapidity. This is 
justified as long as one is restricted to small rapidity volumes \rj\ <C 1. The lattice parameters (all 
of which are dimensionless) are 

• -^Vj.) the number of lattice sites in the transverse direction, 

• Nn, the number of lattice sites in the longitudinal direction, 

• 5 f2 A ta -L) the lattice spacing in the transverse direction, 

• a v , the lattice spacing in the longitudinal direction, 

• Ti n it/a±, the proper time at which the 3-dimensional simulations are begun, 

• St, the time stepping size, 

• A, the initial size of rapidity fluctuations, 

• b, the strength of the "band filter" for the large seed case (see Eq. 1)221) ). 



Of these, only the combinations g 2 fJ,a±Nj_ = g 2 \iL and a^N^ = have a transparent physical 
meaning 5 . Given the energy of the collision and the size of the colliding nuclei, one can estimate 
g 2 [i j^El, since we have periodic boundary conditions, L 2 = ttR 2 . For RHIC collisions of gold 
nuclei, one has g 2 [iL ~ 120; collisions of lead nuclei at LHC energies will be a factor of two larger. 
Also, L v is the physical size of the region in rj being studied. For the purposes of this study, we will 
consider L v = 1.6 units of rapidity in most cases. The continuum limit is approached by keeping 
g 2 fiL and fixed while sending 5t — > 0, g 2 /J>a± — > 0, a v —* 0. For the 3-dimensional simulations, 
we still have to choose a value for Ti n a, which should be such that for A = we stay very close to 
the result from the 2-dimensional simulations (for all of which Ti n u = 0). Thus, we set 

Tinit = 0.05 a ± , (24) 

and later check how strongly the results obtained depend on this choice. 

IV. THE WEIBEL INSTABILITY 

The primary observables injsimulations of classical Yang-Mills dynamics are the components of 
the energy-momentum tensor 



T>» = -g^g uP g^F ai Fp 5 + ^VV 5 *^ (25) 



In particular, 



rpxx , rpyy — 2 Tr 

T 2 T r,r, = T -2 Tr 



F xy + E ?i 



F%i + E i 



Tr 



7?2 I rp2 

xy ' 7] 



(26) 
(27) 



5 The initial size of the rapidity fluctuations A and the band filter parameter b will presumably be further specified 
in a complete theory. For our present purposes, they will be treated as arbitrary parameters, and results presented 
for a large range in their values. 
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FIG. 1: tPl(t, v) as a function of momentum v, averaged over 160 initial conditions on a 16 x 2048 lattice 
with g 2 [iL — 22.6 and L v = 102.4, A ~ 10 -11 . Four different simulation times show how the softest modes 
start growing with an distribution reminiscent of results from Hard- Loop calculations Also indicated 
are the respective values of f max for three values of g 2 /J,r (see text for details). 



Furthermore, the relation H = tT tt = t(T xx + T yy + t 2 T^) (c.f. Eq. ifTSjl ) shows how the physical 
energy density T TT is related to the Hamiltonian density. Similarly, we introduce 

T P ± = T - (T xx + T yy ) , tP l = r 3 T w , (28) 

which correspond to r times the mean transverse and longitudinal pressure, respectively. 

When studying the time evolution of rapidity-fluctuations, it is useful to introduce Fourier 
transforms of observables with respect to the rapidity. For example, 

P~L(T,v,k± = 0) = J drjexp(ir]v){PL(T,x_L,r])}± , (29) 

where denotes averaging over the transverse coordinates (x, y). Apart from v = 0, this quantity 
would be strictly zero in the boost-invariant (A = 0) case, while for non-vanishing A and u, Pl{v) 
has a maximum amplitude for one specific momentum v. Using a very small but finite value of 
A, this maximum amplitude is very much smaller than the corresponding amplitude of a typical 
transverse momentum mode. The system is very anisotropic in momentum space and consequently 
prone to develop a Weibel-type plasma instability ^3], as was shown previously in [25I l3(i|. 

In order to illustrate some features of this instability, it is convenient to excite very low frequency 
rapidity modes. Thus, contrary to our earlier requirement, we have to use very large (Ar] ~ 100) 
rapidity volumes; note that we only use this for illustrative purposes and shall later consider much 
smaller volumes. 

In Fig. ^ we show the ensemble- averaged tPl(v) for four different simulation times. The earliest 
time (g 2 /ur ~ 22) shows the configuration before the instability sets in. Starting to overcome the 
effect of expansion at times q 2 UT — 28, the instability makes the amplitude of the soft modes 
grow like ~ exp(- v /r) (see [23], but also [23] for a more detailed analysis). At times g 2 fir ~ 72 
the amplitudes of the softest modes clearly differs from the starting configuration. The two later 
snapshots (for times g 2 ^r ~ 156 and g 2 \XT ~ 288) indicate that the growth rate of the unstable 



modes closely resembles the analytic prediction from Hard-Loop calculations 13, [la]. Fig. ^ also 
shows that the unstable mode with the highest mode number v (which, in the following we will 
refer to as !/ max ) moves to higher v as a function of time. We will investigate the very interesting 
and suggestive behavior of f max more closely in the next section. 
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FIG. 2: Left: Time evolution of u max , on a lattices with g 2 fJ.L = 22.7, L v — 1.6 and various violations of 
boost-invariance A. The dashed line represents the linear scaling behavior. Right: Time evolution of the 
maximum amplitude tPl (t, v) ; when this amplitude reaches a certain size (denoted by the dashed horizontal 
line), t'max starts to grow fast. 



The unstable mode with the biggest growth rate (the cusp of the "bumps" in FigQ) also seems 
to move upwards in v as a function of r, but more slowly. We use this to define the average growth 
rate T av by simply tracking the time evolution of the maximum amplitude of tPl{v) and fitting it 
with the functional form 

C0 + ciexp(r^ v /^V). (30) 

The label "average" refers to the fact that we do not fit the growth rate of a individual gauge 
mode but rather a convolution of modes contributing to every single v of tPl(v). However, this 
method has the benefit of being a gauge-invariant measurement, while tracking the growth of 
individual gauge modes necessarily involves gauge-fixing procedures which we will address in the 
near future 281. 



A. Results for the growth rate 



Since the components of the energy-momentum tensor are gauge-invariant, the fast rise of the 
soft longitudinal momentum occupation numbers cannot be a gauge artifact. Instead, in principle, 
this fast rise could be due to a lattice artifact rather than a physical instability. To convince 
the reader that this is not the case, we therefore present results in this section for the fitted 
average growth rates obtained for various values of the lattice parameters. The results are shown 
in Appendix FBI and indicate that the fitted average growth rate is close to ~ 0.5, independent of 
our choices for a^, g 2 fia±,L ri , T m it, A and having a weak dependence on g 2 \iL. Our results suggest, 
unequivocally, that the the instability present in the Glasma is genuine and no lattice artifact. 



B. Time Evolution of f max 

Physically, as suggested by Fig. i/ miH denotes the largest mode number that is sensitive to 
the instability. The technical method we use to determine v max , at a given r, is to fit a high order 
polynomial to the momentum spectrum tPl(t, v) and find its section with a pre-defined "baseline" 
(a constant in v and r). Varying the "baseline" by a factor of yfe ~ 1.65 provides us with a rough 
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FIG. 3: Time evolution of the (ensemble-averaged) maximum amplitude of tPl{t,v), for g 2 ^L = 22.6, 

1 /2 

L v = 1.6, iVj_ = 16, 32, A = 0.1 a v ' and N v ranging from 32 to 2048. Larger lattices correspond to smaller 
A. This explains why the early-time behavior is not universal for the simulations shown here. 



estimate of the error of our method. The time evolution of z^max is plotted in Fig [2] for different 
lattices. 

From this figure, one observes an underlying trend indicating a linear increase of ^ max with 
approximately z^ max — 0.06 g 2 [it. For sufficiently small violations of boost-invariance, this seems 
to be fairly independent of the transverse or longitudinal lattice spacing we have tested. 

Presumably, this can be understood as follows: an unstable mode with wavelength Az will cor- 
respond approximately to a wavelength Ary as A77 — atanh-^. At late times (or central rapidities), 
we have t ~ r 3> Az and accordingly 

(31) 

where v ~ (At/) -1 . In other words, the wavelength Ar\ of the unstable mode would decrease (and 
v increase) linearly as a function of r for fixed Az. This is precisely what is seen in Fig. El Turning 
the argument around, Fig. |2] provides us with a measure of the unstable mode wavelength Az for 
large anisotropies, which cannot be calculated within a Hard-Loop framework. 

For much larger violations of boost-invariance - or sufficiently late times - one observes that 
i/ max deviates strongly from this "linear law". In Fig. [2] we show that this deviation seems to occur 
when the maximum amplitude of tPl(t, u) reaches a critical size, independent of other simulation 
parameters. This critical value is denoted by a dashed horizontal line and has the magnitude 
3 • 10~ 5 in the dimensionless units plotted there. A possible explanation for this behavior is that 
the critical size of the longitudinal fluctuations (corresponding to transverse magnetic field modes 
with small fcj_-see Eq. (j27l29|l ) is sufficient to bend "particle" (hard gauge mode) trajectories 
out of the transverse plane into the longitudinal direction. This is essentially what happens in 
electromagnetic plasmas. Note however that in e.m. plasmas, the particle modes are the charged 
fermions, while here the particle modes are hard ultraviolet transverse modes of the field itself. 

In the next section, we will present results that indeed show that the transverse dynamics is 
visibly affected once the longitudinal fluctuations have reached a critical size. However, at the 
present, we cannot exclude the possibility that the rapid rise of z/ max is due to a non-Abelian 
turbulent cascade as described in or the phenomenon described in Clarifying this issue 
will require a detailed analysis of the time dependence of hard and soft modes [2^ . 
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FIG. 4: Time evolution of the ensemble and volume averaged longitudinal pressure Pl, for lattices with 
g 2 \xL = 22.6, L v = 1.6, Nj_ = 16,32, A = 0.1a, 1 / 2 and N v ranging from 32 to 2048. A reduced statistical 
ensemble of only 2 runs for N n = 4096 is consistent with N v — 1024, 2048 results within error bars. 



V. LARGE SEEDS 

In this section we focus on the (more realistic) model that involves large and predominantly 
soft momentum rapidity fluctuations. The only difference with respect to results from the previous 
section is the use of Eq. Q22|) . where - unless otherwise stated - we used the value b = 0.5. 

In FigEl we plot the time evolution of the maximum amplitude of the ensemble averaged 
tPl(t,v), for lattices with different a^. Early times in this figure (g 2 /-iT < 200) correspond to 
the stage when the Weibel instability is operative. Interestingly, all our simulations then show a 
saturation of the growth at approximately the same amplitude. This suggests that what one is 
seeing is similar to the phenomenon of "non-Abelian saturation" jfound in the context of simulations 
of plasma instabilities within the Hard-Loop framework [21I E^ . However, while this result is 
robust for very large longitudinal lattices, further studies involving much larger transverse lattices 
are required to conclusively establish this result. 



A. Creation of longitudinal pressure 



While for the small seed case the longitudinal fluctuations always carried only a tiny fraction 
of the total system energy the situation is different for the large seed case. In the simulations 
reported here, the initial energy contained in the longitudinal fluctuations is about 1% of the total 
system energy. In this subsection, we study the behavior of the time evolution of the average Pl(j]). 
This value is consistent with zero in simulations that do not allow for longitudinal dynamics. In 
Fig. ^ we plot Pl as a function of r for different lattice spacings a„. As can be seen, for large 

(meaning low lattice UV cutoff), the longitudinal pressure is consistent with zero while it is 
clearly non-vanishing when the lattice UV cutoff is raised. However, there seems to be a limit to 
this rise as there is no notable difference between the simulations for the three smallest values of 
the lattice spacing. This is suggestive that the rise in the longitudinal pressure is physical, though 
again, further studies on even larger transverse lattices are needed to strengthen this claim. 
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FIG. 5: Volume and ensemble averaged Hamiltonian density {H{r,rf)) n (corresponding to r times energy 
density) and transverse and longitudinal pressure t{Pt(t, r))} 71 , t(Pl(t, for simulations with small 
(16 x 32 lattice) and large (16 x 2048 lattice) longitudinal UV cut-off; quantities are multiplied by t to 
facilitate comparison. The system evolves towards a more isotropic state if the UV cut-off is sufficiently 
large. A fit shows that the energy density behaves as e ~ T - 1 067 a t late times. All curves calculated on 
lattices with g 2 [iL — 22.6, L v = 1.6 and A = O.la^ 2 . 



B. Towards isotropy 

In Fig. |5J we investigate the time evolution of the transverse pressure as well as the energy density 
for (i) a simulation with a low UV cutoff (16 x 32 lattice) and (ii) a simulation with a high UV 
cutoff (16 x 2048 lattice). We observe that the rise in the mean longitudinal pressure accompanies 
a drop both in the mean transverse pressure and energy density, thus bringing the system closer to 
an isotropic state. The energy density depends on the proper time as £ ~ i {^ 7 , which is distinct 
from e ~ for an isotropic system. Despite this clear trend, no full isotropization was achieved 
in these simulations on the time scales under investigation. 

Why is this so? We digress here to emphasize, as discussed previously, that while simulations of 
the classical equations of motion capture important aspects of the early time dynamics, they miss 
others. What contributions are included in the numerical simulations here, and what contributions 
are not, can be understood in a systematic formalism developed recently to treat the real time 
non-equilibrium dynamics of fields in the presence of strong (j ~ |) time dependent sources (2r3 |. 
The power counting in these theories, as discussed in Ref. [2f|, is entirely in powers of g 2 h. In 
this power counting, the 2+1-dimensional boost invariant classical solution provides the leading 
contribution of order l/g 2 h to the average particle multiplicity- but includes all orders in the 
sources (g j ~ 0(1)). This leadin g or der result was computed numerically by Krasnitz, Nara and 
Venugopalan @, 0| and by Lappi [Tnl] . At next-to-leading order (NLO), which is of order (g 2 h)°, 
there are two contributions to the average multiplicity: I) from the small fluctuations propagator, 
and II) from the product of the classical field at lowest order and the classical field at one loop 
order. NLO contributions of type I are obtained by solving partial differential equations (with 
retarded boundary conditions) for small fluctuations on top of the classical background field. The 
procedure followed here is very similar in spirit: we add fluctuations 5Ei and SE^ on top of the 
boost invariant classical solutions and study their evolution in time. Whether these are identical to 
the NLO contributions is a topic under active investigation and beyond the scope of this paper (but 
which will be addressed further in a forthcoming work [13|). However, it is clear that contributions 
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similar to the type II NLO contributions are not included in our simulation. 

It has been argued ji^] that the early time classical dynamics can be smoothly matched on to 
a Boltzmann equation. Interestingly, the previously mentioned NLO contributions may be clearly 
identified in the Boltzmann equation correspondence of Ref. - such a study would shed further 
light on whether the additional processes not included in the study here may play a role in speeding 
up thermalization. 

In any case, since the present study is based on classical field dynamics only, reaching true 
thermal equilibrium is beyond the capabilities of our simulation since description of the late stages 
of the evolution necessarily require the inclusion of quantum effects. Nevertheless, we believe that 
our study is at least qualitatively valid until the onset of the equilibration process, e.g. the point 
where the energy density has departed from the free-streaming behavior e ~ r _1 . As argued above, 
adding quantum effects through the full inclusion of the NLO contributions may be used to extend 
the domain of applicability of such simulations. 

We return to our discussion to note, however, that increasing the seed parameter A further 
pushes the trend toward isotropization to earlier times. In Fig. H3 we plot the time evolution 
Pt/Pl as a measure of the system anisotropy for three values of A. One observes that larger 
values of A accompany an earlier set-in of the isotropization process. It should also be noted that 
increasing A further has the effect of strongly increasing the total system energy at early times. 

Similar to increasing A, one can also decrease the value of b (where we recall that this is the 
parameter that controls the initial exponential fall-off of the momentum spectrum). Decreasing b 

1 /2 

by a factor of 5 results in a system anisotropy that is similar to the values for A = 10 a v in Fig. |BJ 
We conclude here that a first principles analytic computation of A and b would be extremely 
useful because quantitative results for the system isotropization time scale appear to depend on 
them strongly. 



VI. CONCLUSIONS AND OUTLOOK 

We discussed here detailed results from numerical simulations of 3+1-D SU(2) Yang-Mills equa- 
tions for an unstable Glasma expanding into the vacuum after a high energy heavy ion collision. 
This work greatly expands on our earlier letter j2j| where first results for 3+1-D numerical sim- 
ulations were presented. Specifically, we explained in detail the set of initial conditions we are 
studying, providing details of our numerical simulations, and tested that our results were robust 
in the longitudinal and transverse continuum limits. A novel addition to our previous studies is a 
detailed study of large violations of boost invariance. 

We established that for Color Glass Condensate initial conditions, a Weibel instability is present 
in the system that appears as rapidly growing longitudinal fluctuations (starting in the far infrared) 
from times of g 2 ^T ~ 30 onwards. We performed detailed checks on the largest available lattices 
that the measured growth rate seems to be nearly independent of all our lattice parameters. We 
have demonstrated that the presence of the instability is very unlikely to be a lattice artifact. 

We investigated the distribution of the unstable modes and found it to be strikingly similar to 
the analytic results from Hard-Loop calculations 0,1110. The time evolution of the maximum 
amplitude of these modes indicate a saturation of the growth at a certain size. This effect appears 
to be identical to the phenomenon of non-Abelian saturation described within the Hard-Loop 



framework |21|, |22j. Tracking the time evolution of the hardest unstable mode, f max , we deduced 
that the smallest unstable wavelength Az is finite even for extreme anisotropies. This goes beyond 
results from Hard-Loop techniques. 

We find further that when the longitudinal fluctuations reach a critical size, f max increases 
dramatically filling up the mode spectrum. Because this effect is simultaneously accompanied 
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FIG. 6: The mean system anisotropy as measured by Pt/Pl for various values of A. Note that the late 
times values could be subject to finite lattice spacing artifacts. 

by a drop in the transverse pressure, we interpret this result as the effect of the Lorentz force 
exerted by modes of the transverse magnetic fields on the hard transverse gauge modes. These 
fields effectively bend the hard transverse gauge modes into the longitudinal direction. Note that 
the total amplitude of the transverse magnetic field does not increase significantly. Nevertheless, 
the "bending" can be accomplished since the instability has populated predominantly modes with 
k± = and the corresponding transverse magnetic fields thus can act over long transverse distances. 

Simulations on lattices with large longitudinal UV cutoff suggest that significant longitudinal 
pressure is built up in this process. This therefore shows a clear trend towards isotropization of 
the system. It is also reflected in the time evolution of the energy density. The energy density 
changes from the free streaming e ~ r _1 behavior to e ~ r" 1 067 . Our results are therefore proof in 
principle that a trend towards isotropization within a weak coupling framework is indeed possible, 
contrary to recent claims 41]. 

However, the time scales associated with isotropization in our simulations, are much too large 
in order to be of interest for phenomenology in heavy-ion collisions. This may be due to the fact, 
as discussed in section V.B, that we are only partly including next-to-leading order corrections in 
our framework [2^| . While these contributions are formally NLO and may be presumed small, they 
properly treat number changing processes. In this regard, they may be considered leading order 
effects in driving the system towards equilibration. We should note further that our results strongly 
indicate that isotropization time scales are much shorter if longitudinal fluctuations are already 
large initially than if they have to be built up by the instability. A first principles calculation of the 
longitudinal gluon spectrum after a heavy-ion collision would thus probably be a key ingredient in 
understanding the process of equilibration. 
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APPENDIX A: LATTICE DISCRETIZATION 

Introducing the lattice gauge links Ui = ex.p(iga±Ai), U v = exp^ga^A^) , and making the fields 
dimensionless by scaling out the transverse lattice spacing a± and the coupling for convenience, 

i*-— , Ai — ► — > ^-4- ( A1 ) 

the discretized form of the Hamiltonian density Eq. (|12|) becomes the analog of the Kogut-Susskind 
Hamiltonian, 



7~i L = — ^-Tr 
rg l a\ 



S i □ 



(A2) 



where we also scaled the lattice time r — ► tl = r/a±; Un = Uij with the standard plaquette 
U^ v {x) = U fJt (x)U l/ (x + [i)UUx + L^J (a?) . The lattice equations of motion then become 

E?(tl + ^,x) = Ei{r L - ^L,x) + 2i5r L t l Tr ( T a Ui(r L , x) E Sj^x) 



where 



2 " ' /I- 

+2i^Tv [r a Ui{T L ,x) E ^(r L ,x 
TLa ri V 

+ = E v (r L - S -^,x) + 2i-^Tr\r a U v (T L ,x) E S^fa, 

Ui{T L + 5t l ,x) = exp(i5T L t^ 1 Ei(r L + ^,x))Ui(T L ,x) 

Stl 

U v (t l + 5t l ,x) = ex.Y>(i5T L T L a rj E ri (T L + —,x))U v (t l ,x) (A3) 



Sl v {r L , x) = U v {t l ,x + (i)Ul(r L ,x + ^[^(tl, z) (A4) 



is the gauge link staple. For SU(2), E % = E % a r a with r a = ^a a and a a the Pauli matrices and 
Tr(r a T b ) = ^<5 afc . Note that the sum Su-i^j runs over both positive and negative directions, with 

U.j{T L ,x) = U](x-j). 

Numerical stability requires a very small 5tl at the beginning of the simulation, but not at late 
times. Thus, we have found it convenient to use adaptive time steps, choosing 

Sr L = e-Z^, (A5) 

with T ~ 20 and e ~ 0.025 giving satisfactory performance for most simulations. 

1. Initial conditions on the lattice 

For SU(2), the initial conditions Eq. (|18jl lead to the following form of the gauge links 

Ui = (^ (1) +^ 2) )K )t+ ^ (2) T 1 ( A6 > 
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U v = 1 2 (A7) 

u\ lU2 \x) = vW>W(x)vW>W(x + i) (A8) 

vW\x) = exp(iAi )2 (a;)) (A9) 

A L Ai, 2 (x) = -Pi, 2 (x) (A10) 

< pUx)fi(y) > = 9 4 p 2 al5 ab 6l y (All) 

<~pl{x)p b 2 {y)> = 0, (A12) 



where x, y here denote site indices in the transverse plane and there is no summation over i on 
the l.h.s. of the first equation. The last two equations should be understood such that the gauge 
configurations pi j2 are drawn as Gaussian random numbers with weight g 2 pa± while there should 
be no correlation between the different nuclei. Finally, Al is the lattice Laplacian in the transverse 
plane, 

A L A(x) = -2A(x) + A (x + i) + A(x - i). (A13) 

i 

The initial condition for the momenta for the boost-invariant case are 

E% = *Tr r a J2 (v}%) + U^(x)uj(x) 

i 

-U\ 2 \x -i)-U}{x- i)U^\x - ijj - h.c. (A14) 
together with Ef(x) = 0. The rapidity fluctuations are constructed as 

5E; = -F(x v )DfSEl(x ± ) 

5E?{x) = (F(x ri )-F(x ri -Ti))6E?(x ± ) 

< 5E?( Xi _)5E b {y L ) > = SPft^Jij 

< F(x v )F(x' v ) > = a-'AX^- (A15) 

For the case of large seeds, we modify F(x v ) by creating its Fourier-transform w.r.t. rapidity, 

F(k) = ^2exp(2-iriKx v /N v )F(x v ), (A16) 

removing the high-frequency modes and transforming back, 

F(x v ) = j-r ^2 exp(-2mKx v /N v )exp(-2Trb\K\/N v /a TI )F(K) 

iVr ' n<N n /2 

+-jrr X! ex P( — 2mKx v /N v ) exp(27rfe(|ft| - N v )/N v /a v )F(n). 

Nr > K>N n /2 

APPENDIX B: MEASURED GROWTH RATES 



We present here a collection of tables with measured growth rates for various different lattice 
parameters. 



N v 


"nan 
1 fit 


32 


0.50 ± 0.01 


64 


0.52 ±0.02 


128 


0.52 ±0.02 


256 


0.524 ±0.016 



TABLE I: Testing for the rapidity lattice 
spacing dependence on lattices with g 2 \iL = 
22.6,L V = 1.6 with N ± = 16 and A = 

10- 10 af. 
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N v 


T^av 

1 at 


32 


0.5 ± 0.01 


64 


0.508 ± 0.012 


128 


0.488 ± 0.01 


256 


0.492 ± 0.016 


512 


0.543 ± 0.016 



TABLE II: Testing for the rapidity vol- 
ume dependence on lattices with g 2 fiL = 
22.6, a v = 0.05 with N ± = 16 and A = 

io- io a y 2 . 



N ± 


N v 


-pav 
1 fit 




N ± 


N v 


g 2 ^L 


-pav 
1 fit 


16 


32 


0.50 ±0.01 




16 


32 


22.6 


0.50 ±0.01 


32 


64 


0.504 ± 0.01 




64 


64 


67.9 


0.426 ±0.01 


64 


64 


0.488 ± 0.02 




128 


64 


90.5 


0.45 ± 0.03 










128 


256 


181 


0.394 ±0.04 



TABLE III: Testing for the transverse lattice 
spacing dependence on lattices with g 2 fiL = 



22.6,L„ = 1.6 and A = 10 



-10 1/2 



TABLE IV: Testing for the transverse vol- 
ume dependence on lattices with L„ = 1.6 



and A = 10 



-10 1/2 
U77 



Tinit/a_L 


-pav 
1 fit 


0.025 
0.05 
0.1 


0.488 ± 0.014 
0.504 ±0.01 
0.51 ±0.01 



TABLE V: Testing for the dependence on 
Tinit on lattices with g 2 fiL = 22.7, L v = 1.6, 
N± = 32, N v = 64 and A = 10- w a v /2 . 



-log 10 (Aa^ 1/2 ) 


-pav 
1 fit 


10 


0.504 ±0.01 


6 


0.514 ±0.03 


3 


0.48 ±0.05 



TABLE VI: Testing for the dependence on 
the initial seed amplitude A on lattices with 
g 2 fiL = 22.7, L v = 1.6, N_l = 32, N v = 64. 



e 


-pav 
1 fit 


0.025 
0.0125 
0.00625 


0.50 ±0.01 
0.50 ±0.01 
0.51 ±0.01 



TABLE VII: Testing for the dependence on the time step St lattices with g 2 [iL = 22.7, L v = 1.6, 
N± = 16, N v = 32. Note that adaptive step sizes are used (see Eq. (|A5[0 . 
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